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ABSTRACT 

The surface-brightness profiles of galaxies are well described by the Sersic law: systems with high 
Sersic index m have steep central profiles and shallow outer profiles, while systems with low m have 
shallow central profiles and steep outer profiles. R. Cen (2014, ApJL, 790, L24) has conjectured 
that these profiles arise naturally in the standard cosmological model with initial density fluctuations 
represented by a Gaussian random field (GRE). We explore and confirm this hypothesis with Wbody 
simulations of dissipationless collapses in which the initial conditions are generated from GREs with 
different power spectra. The numerical results show that GREs with more power on small scales lead 
to systems with higher m. In our purely dissipationless simulations the Sersic index is in the range 
2 < m < 6.5. It follows that systems with Sersic index as low as m ~ 2 can be produced by coherent 
dissipationless collapse, while high-m systems can be obtained if the assembly history is characterized 
by several mergers. As expected, dissipative processes appear to be required to obtain exponential 
profiles {m ~ 1). 

Subject headings: galaxies: bulges — galaxies: elliptical and lenticular, cD — galaxies: formation — 
galaxies: fundamental parameters — galaxies: structure 


1. INTRODUCTION 


The surface-brightness profiles of spheroidal an d disk 
compo nents of galaxies are well described by the iSersid 

(CSlaw 


I (R) = Iq exp 



( 1 ) 


where Re is the effective radius, m is the Sersic index, 
le = I {Re) is the surface brightness at the effective ra¬ 
dius and b{m) 2m — 1/3 + 4/(405m) (|Ciotti fc BertinI 
I1999D . When m is high the central profile is steep and the 
outer profile is shallow, while when m is low the central 
profi l e is sh allow and the outer profile is steep. Recently 
iCenI (|2Q14l ) has envisaged that these profiles arise nat¬ 
urally in the standard cosmological model with initial 
density fluctuations represented by a Gaussian random 
field (GRE). The underlying idea is that central concen¬ 
trations of stars and extended envelopes are formed by 
the late infall and accretion of substructures. Therefore, 
if the fluctuation field is dominated by long-wavelength 
modes the formation of the galaxy, due to the absence 
of significant substructures, is mainly determined by a 
coherent collapse, and the final profile will be shallow in 
the center and steep in the outskirts (an extreme case is 
the exponential profile m = 1, which is typical for disks). 
If the fluctuation field is dominated by short-wavelength 
modes there is substantial late infall of substructures and 
the final profile will be steeply rising toward the center 
and gently declining in the outskirts (as observed, for 
instance, in massive elliptical galaxies with m > 4). 

In this Letter we explore quantitatively Gen’s pro¬ 
posal with numerical experiments in which we follow 
the dissipationless collapse of cold distributions of par¬ 
ticles whose initial conditions are determined by GREs 
with different power spectra. Numerical simulations of 
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dissipa tionless collapse have been run by several au¬ 
thors (Ivan Albad^ 119821: lUdryl 11993: iBoiiv et al.l 2002 


__ __ __ 

Nipoti et al.l 120061: 1 Joyce et al.l 120091: iSylos-Labinil 2013 : 

Benhaiem fc Svlos Labinill20l5l : 1Worrakitpoonponll20l5[) . 


A general finding is that the end-products have surface- 
density profiles well fitted by the Sersic law (equation [T]). 
Remarkably, this result is not specific to Newtonian grav¬ 
ity, as it is also found in studies of dissipationless col¬ 
lapses in modified gr avity theories (jNipoti et al.l I2007I: 
iDi Gintio et al.|[2013f ). In Newtonian gravity the final 


Sersic index is typically c lose to m 


lAguilar fc Merrit 


- 

responding to the lde VaucouleursI (|l948l ) profile. Though 


hereafter 


4 (Ivan A]ha,flal[T^ 


AM90f) . the value cor- 


there are some indications that the final value of m 
can depend on the initial conditions (jTrenti et al.ll2QQ^ 
iNipoti et al.ll2QQ^ . so far there is no clear evidence of 
a dependence of m on the propertie s of the initi a l fiuc - 
tnations. Since the seminal work of Ivan Albadal (ITo^ 
it was realized that the dumpiness of the initial condi¬ 
tions is an important factor in determining the nature of 


considered in several investigations ( 

Mav & van Albada 

1984 

: iMcGlvnnl11984 iLondrillo et al. 

119911: iRpv & Perez 


homogeneities of the initial phase-space distribution 
were not systematically classified in terms of fl uctuation 
powe r spectrum. An exception is iKatzl (|1991l . hereafter 
lK91f ). who set up the initial conditions self-consistently 
fro m a GRE and explored diffe rent power spectra (see 
also iDubinski fc Cariber^ll991l and iBinnev fc Tremaine! 
I2QQ8I. section 4.10.3b 

Here we present high-resolution A^-body simulations 
aimed at isolating the effect of the spectrum of the inho¬ 
mogeneities of the initial conditions on the final density 
profiles of cold dissipationless collapses. We present ev¬ 
idence that the Sersic index of the collapse end-product 
correlates with the slope of the fluctuation power spec¬ 
trum of the initial conditions. 
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TABLE 1 

Properties of the simulations. Pq and n: amplitude and 

POWER-LAW INDEX OF THE FLUCTUATION POWER-SPECTRUM OF 
THE INITIAL CONDITIONS . c/a AND b/a: FINAL 
SHORTEST-TO-LONGEST AND INTERMEDIATE-TO-LONGEST AXIS 
RATIOS. Thalf/'^O: FINAL HALF-MASS RADIUS IN UNITS OF THE 
SCALE RADIUS Tq. m AND am- BEST-FITTING SeRSIC INDEX OF THE 
FINAL DENSITY PROFILE AND ASSOCIATED UNCERTAINTY. 


Name 

Po 

n 

cja 

bja 

^half/^O 

m 

m 

PO 

0 

_ 

0.43 

0.56 

0.92 

2.01 

0.04 

P03n3 

0.3 

-3 

0.46 

0.71 

1.05 

2.26 

0.02 

P03n25 

0.3 

-2.5 

0.46 

0.76 

1.04 

2.23 

0.02 

P03n2 

0.3 

-2 

0.48 

0.83 

1.08 

2.38 

0.02 

P03nI5 

0.3 

-1.5 

0.53 

0.91 

1.05 

2.61 

0.04 

P03nl 

0.3 

-I 

0.61 

0.98 

1.01 

3.32 

0.12 

P03n05 

0.3 

-0.5 

0.60 

0.93 

0.68 

4.90 

0.15 

P03n0 

0.3 

0 

0.63 

0.92 

0.50 

6.43 

O.II 


2. NUMERICAL EXPERIMENTS 
2.1. Initial conditions 

The initial conditions of each simulation are built as 
follows. Working in Cartesian coordinates, we take a 
cube of edge /q = 2ro centered in the origin, in which we 
generate a GRF 

<5(x) = ^ E (2) 

k 

where and (5_k are independent random variables with 
(5_k = <^k 5 ^ = (n and the sum is perforni ed over half of 
the k-space (see lBinnev fc Tremainel[2QQM section 9.1.1). 
The GRF S{x.) is fully characterized by its power spec¬ 
trum P(k) = (|(^kP)/R, which we parameterize as 

^ if ko<k<kcntoi^ 

^ ^ " lo if k<ko or k> A^cutoff. ^ ^ 

Here Pq and n are, respectively, the power-spectrum am¬ 
plitude and index, ko = 27vjro is the minimum wave- 
number and /^cutoff is the cut-off wave-number. 

We build a lattice of Nx x Ny X Nz points equally 
spaced in x, y and z. At each lattice point we compute 
the density fielcfl p(t, z) = pbg exp((5), where Pbg(^) = 
po{r/ro)~^ is an unperturbed background power-law dis¬ 
tribution (r = and 5{x^y^z) is the GRF 

defined above. Finally, p is normalized to its maximum 
value, so 0 < p < 1 independent of po- 

Using the standard rejection technique, we generate 
a spherical distribution of N equal-mass particles with 
r < ro and density distribution p. The t, y and z coor¬ 
dinates of each particle are then multiplied, respectively, 
by Qy and Qz to get a triaxial configuration. Velocity 
components Vx, Vy and Vz, extracted from a Gaussian 
distribution with vanishing mean and unit variance, are 
temporarily assigned to each particle. Once the total po¬ 
tential energy W and the temporary total kinetic energy 
T are computed, all the velocity components are mul¬ 
tiplied by ^/^\W\p2T in order to obtain a system with 
initial virial ratio jS. 

^ We assume that the amplitudes of the fluctuations generated 
by the GRF are distributed log-normally because our initial con¬ 
ditions are meant to represent the non-linear phase of the collapse 
(IKavo et al.l[200TlL 


2.2. Parameters of the simulations 

The aim of our simulations is to isolate the effect on 
the end-product of dissipationless collapse of the relative 
contribution of short- and long-wavelength modes in the 
fluctuation power spectrum. Therefore, we present the 
results of a set of simulations differing only in the value of 
the power-spectrum index n, which spans the range — 3 < 
n < 0 (see Table [T]). In all the simulations the phases 
of the GRF modes are the same, as well as the other 
parameters determining the initial conditions: Pq = 0.3, 
^cutoff = 6/co, 7 = 1, /3 = 0.01, Qx = 6, Qy = 4, qz = 2, 
and Nx = Ny = Nz = SO. For comparison, we also ran 
a simulation (named PO) with the same parameters as 
above, but no fluctuations (Pq = 0). 

The V-body simulations were run with the parallel col- 
lisionl ess code fvfps (fortran version of a fa st Poisson 
solver: iLondrillo et al.ll2003l: iNipoti et al.ll200^ . which is 
based on iDehnei] (|2002f ) algorith m, an efficient combina¬ 
tion of the fast multiple method (iGreengard fc R^okhlinl 
119871 ) and the tree code (jBarnes fc Hntlll986[ ). The main 
parameters of the fvfps code are the number of parti¬ 
cles V, the softening length e, below which the Newto¬ 
nian force is smoothed, and the minimum value of the 
opening parameter ^min = 0.5, which determines the 
mass-dependent tolerance parameter 0 (analogous to the 
opening angle of iBarnes fc Hntiri986[ ) used to control the 
acc uracy of the force ap proximation (see iDehnenI l2QQ2l 
and lLondrillo et aD2QQ3[ ). In our simulations we adopted 
N = 10^, 6>min = 0.5 and 5 = 0.02ro. The time-step At, 
which is the same for all particles, is allowed to vary 
adaptively in time as At = 0.3/(47rGpmax)^^^, where 
Pmax is the maximum particle density. Each simulation 
runs from t = 0 to t = lOOtu, where tu = a/fq /GM 
and M is the total mass of the V-body system. With 
this choice the system is fully virialized at the end of the 
simulation (the initial system’s free-fall time is « lOtu). 
The time-step values are in the range 8 x 10-4 < At/tu < 
3 X 10-2. 

2.3. Results 

The intrinsic and projected prop erties of the c ollaps e 
end-products are determined as in INipoti et al.l (j2QQ6f ). 
At the end of the simulation only bound particles are 
selected (the mass loss is in the range 2-11%). The po¬ 
sition of the center of the system is determined usin g 
the iterat ive technique descri bed by iPower et al.l (|2QQ3f) . 
Following INipoti et al.l (|2QQ2[) . we measure the axis ratios 
c/a and 6/a of the inertia ellipsoid of the final density dis¬ 
tribution, its angle-averaged profile and half-mass radius 
Fhaif {cl, b and c are, respectively, the longest, intermedi¬ 
ate and shortest axes). 

In Fig. [1] we plot the initial and final particle distribu¬ 
tions (projected along the shortest axis), together with 
the final angle-averaged density profiles of four repre¬ 
sentative simulations: simulation PO (with Pq = 0) and 
simulations P03n3, P03nl and P03n0 (with Pq = 0.3 
and, respectively, n = —3, n = —1 and n = 0). As ap¬ 
parent from the plots in the left-hand column of panels, 
the initial distribution becomes more and more clumpy 
from top (simulation PO, with homogeneous initial condi¬ 
tions) to bottom (simulation P03n0, in which most par¬ 
ticles are initially in relatively few small clumps). The 
corresponding end-products (central column of panels of 
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Fig. 1. — Initial (left-hand column of panels) and final (central column of panels) particle distributions (projected along the shortest 
axis) of simulations PO, P03n3, P03nl and P03n0 (from top to bottom). The panels in the right-hand column show the corresponding final 
angle-averaged density profiles (circles) together with their best-fitting deprojected Sersic profiles (solid curves); using the right-hand axes 
the circles can be interpreted as dynamical time tdyn as a function of radius for the final systems. Here times are in units of tu, x and y 
are in units of the scale radius ro, and Phalf is the density at the half-mass radius r^aif- 
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Power-spectrum index n 


Fig. [U tend to be more extended when the initial con¬ 
ditions are more homogeneous (upper plots) and more 
compact when the initial conditions are clumpier (lower 
plots). The final half-mass radius rhaif ranges from vq 
when n = —3 down to 0.5ro when n = 0 (see Tabled] 
and Fig. [2^). The end-products are typically triaxial 
with 0.43 < c/a < 0.63 and 0.56 < b/a < 0.98. The 
trend with n is that the systems tend to be almost pro¬ 
late for smooth initial conditions and almost oblate for 
clumpy initial conditions (see Tabled) and Fig. [2]3). 

The panels in the right-hand column of Fig. [T] 
show the final angle-averaged density profiles of the 
aforementioned representative simulations for —1.5 < 
log(r/rhaif) < 1- Over this radial range the dynami¬ 
cal tim43 ^dyn(^) is shorter than the simulation time- 
span lOOtu (see rightmost axes in Fig. [T]), so the den¬ 
sity profiles can be considered stationary. From t hese 
plots it is apparent that, consistent with the idea of iCenI 
(|2Q14l ). the final angle-averaged density profile is steeper 
in the center and shallower in the outskirts when the ini¬ 
tial conditions are clumpier. We quantify this finding 
by comparing the final distributions with the Sersic law 
(equation d]) • Under the assumption of spherical sym¬ 
metry and position-independent mass-to-light ratio, the 
intrinsic stellar mass density distribution corresp onding 
to eq uation ([T]) can be obtained in integral form (jCiotti 
119911 ). A simple approximation of the deprojected Sersic 
profile, which we adopt in this work, is 


p(r) — Phalf 





(4) 


(jLima Neto et al.1 Il999f ). where u = 1/m, p = 1 — 
0.6097z/ + 0.05463z/^, phaif = p(^haif) is the density at 
the half-mass radius and rg is a characteristic radius re¬ 
lated to rhaif by 




Power-spectrum index n 


Fig. 2.— Final half-mass radius rhaif normalized to the scale ra¬ 
dius ro (panel a), axis ratios (panel b) and best-fitting Sersic index 
m (panel c) as functions of the initial fluctuation power-spectrum 
index n of the A/^-body simulations P03n3, P03n25, P03n2, P03nl5, 
P03nl, P03n05 and P03n0 (from left to right; see Tabled}. 


The final angle-averaged density profiles of our Wbody 
simulations are very well represented by the deprojected 
Sersic law (equation S]) . The profiles are fitted taking m 
as only free parameter, because rhaif and phaif are fixed 
by the measured values. The fits performed over the ra¬ 
dial range 0.04 < r/rhaif <10 give values of the Sersic 
index in the interval 2 < m < 6.5 with small associated 
uncertainties 0.02 < am ^ 0.15 (see Table[T]and panels in 
the right-hand column of Fig. [1]). In simulation PO (with 
smooth initial conditions) we perform the fit over the 
smaller radial range 0.04 < r/rhaif < 5, because the pro¬ 
file has a power-law tail at large radii (see top-right panel 
of Fig. [T]), which is reminiscent of the core-halo structure, 
a well-known feature of collapses starting from homoge¬ 
neous initial conditions (jLondrillo et al.lTl991f ). Fig. [2t 
shows that the best-fitting Sersic index m increases for 
increasing n: clumpier initial conditions lead to higher 
values of m. 

^ At each angle-averaged radius r we define the dynamical time 
as tdyn('^) = -\/37r/16Gp(r), where p{r) = 3M(r)/47rr^ is the aver¬ 
age density within r and Mir) is the mass contained within r. 
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3. DISCUSSION AND CONCLUSIONS 

The r esults of our simulations confirm the conjecture 
of iCenI (|2Q14l ): the density profile of dissipationless col¬ 
lapse is steeper in the center and shallower in the outer 
parts if the fluctuation power spectrum of the initial con¬ 
ditions is dominated by short-wavelength modes. Vice 
versa, power spectra dominated by long-wavelength fluc¬ 
tuations lead to density profiles that are shallow in the 
center and steep in the outskirts. The end-products 
of our simulations have density distributions well rep¬ 
resented by the deprojected Sersic law with index in the 
range 2 < m < 6.5. For increasing spectral index n the 
best-fitting Sersic index m increases, the half-mass ra¬ 
dius Thaif decreases, and the systems tend to move from 
almost prolate to almost oblate intrinsic shape. 

Of course, the exact values of the measured quanti¬ 
ties are expected to depend on the details of the initial 
conditions: for instance, while here we find m 2 for 
the end-product of simulation PO (with smooth initial 
p (X r~^ density profile), it is well known tha t the end- 
produ ct of a cold collapse with smooth initial iPlummerl 
(1191 Ilf density distributio n is extremel y well fitted by 
theE e Vaucouleurd (Il948lf m = 4 profile (jLondrillo et al.l 
119911: iNipoti et al.l l2QQ6f ). Therefore, the above range 
2 < m < 6.5 must not be taken at face value. However, 
it is interesting to notice that, based on the results of the 
present work and of previous studies, it appears hard to 
get m < 2 with purely dissipationless processes, consis¬ 
tent with the expectation that the formation of m ~ 1 
systems (typically disks) requires dissipative processes. 

It is interesting to compare our results with those of 
previous simila r inv estigations. A very interesting work 
is the paper of lK91l . who attempted a systematic study 
of the effect of GRF power spectrum on the structu re of 
virialized systems. Though the initial conditions of lK91l 
simulations, which were meant to represent conditions 
before turn-around, are different from ours in many re¬ 
spects, based on the results of the presen t work we should 
expect that the final density profiles of lK91l depend on 
the initia l fluc tuation power spectrum. In fact, the con¬ 
clusion of lK91l is that the final profiles do not depend sig¬ 
nificantly on the power-spectrum slope: however, when 
compared t o tod ay’s standard, the resolution of the sim¬ 
ulations of lK91l is rather poor 4000 particles), so it 
is likely that detailed differences in the density profiles 
were obscured by numerical noise. 

A set-up in a sense more similar to ours was that of 


IAM90I . who did not included fluctuations in their initial 
conditions, but considered the dissipationless collapse of 
smooth triaxial particle distribution with the same initial 
density field as our background distribution (pbg oc r“^). 
Our simulation PO (with no flu ctuation ; Pq = 0) is 
therefore very similar to those of lAM9(j| with virial ra¬ 
tio (3 ~ 10“^ and actually, consistent IaM 90I . we find a 
prolate final system with axis ratios c/a ^ h/a ~ 1/2. 
Howev er, while we find best-fitting Sersic index m 2, 
IAM90I report that their final distributions are well fitted 
by m = 4. Again, this is likely a matter of resolution: 
IAM90I . with 5000 particles, can follow the profile down 
to ~ 0.5rhaif, where the difference between the m = 2 
and m = 4 profiles is hard to detect, while in this work 
we have been able to fit the profiles down to 0.04rhaif- 

Independent support to Gen’s model and to our results 
comes from numerical studies of dissipationless galaxy 
mergers. A galaxy growing in a region of the Universe 
dominated by fluctuations on small scales is expected to 
form by several mergers of smaller subunits. The finding 
that the best-fitting Sersic index m increases for increas¬ 
ing fluctuation power-spectrum index n is therefore con¬ 
sistent with the results of numerical simulations showing 
that d issipationless merg ers make the Sersic index in¬ 
crease (|Nipoti et al.ll2QQ^ . In particular, the dissipation¬ 
less accretion of small satellites is believed to b e the most 
promising mechanism to form high-m systems (|Hilz et al.l 
I2Q13I ). However, dissipative processes can also contribute 
to raise m, as found for instance in s imulations of mergep 
between gas-rich disk galaxies (e.g. lHopkins et aLll2QQ9[ ). 

In this work we have provided quantitative support to 
the idea that the origin of the Sersic law is related to 
the fluctuation power spectrum in the initial conditions 
of galaxy formation. Still, our results are based on toy 
models that neglect all the complexities of proper galaxy 
formation theories. In the future it will be interesting to 
explore this idea more realistically by using cosmologi¬ 
cal simulations with distinct baryonic and dark matter 
components, and including the all-important dissipative 
processes. 
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